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SUMMARY 

v 

The  author's  method  for  computing  steady,  inviscid,  subcritical  flow  past 
a thick  cambered  wing  is  extended  to  design  applications.  Four  problems  are 
considered:  (1)  given  thickness  and  doublet  (first-order  loading)  distributions; 
(2)  given  thickness  and  upper-surface  pressure  distributions;  (3)  given  loading 
and  upper-surface  pressure  distributions;  (4)  a hybrid  of  (2)  and  (3)  in  which 
the  thickness  is  specified  everywhere  except  near  the  root,  and  is  determined 
near  the  root  when  the  doublet  distribution  is  constrained  to  exhibit  spanwise 
invariance  in  that  region.  Convergence  for  the  first  problem  is  excellent. 

For  all  problems,  good  convergence  is  obtained  outboard.  For  the  single  case 
reported  of  the  second  problem,  convergence  was  secured  near  the  root  but  cannot 
yet  be  guaranteed.  Near  the  root,  slow  convergence  was  obtained  for  the  third 
problem,  rather  better  convergence  for  the  fourth  problem.  This  hybrid  option 
is  tentatively  recommended. 


* Replaces  RAE  Technical  Report  76027  - ARC  36857. 
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INTRODUCTION 


This  Report  describes  the  development  of  computer  programs  to  design 

camber  and  twist  distributions,  and  in  some  applications  thickness  distributions 

also,  for  a wing  without  dihedral  in  steady  inviscid  incompressible  (or  sub- 

critical)  flow.  Tht?  programs  are  extensions  of  the  author's  direct-calculation 

program,  described  in  Ref.l,  in  which  the  flow  field  is  represented  by  suitable 

source  and  doublet  distributions  in  the  wing  chordal  surface.  A computer  sub- 

2 3 

routine  is  available  * to  calculate  the  corresponding  flow  fields.  Linear 
compressibility  effects  are  accounted  for  by  working  in  the  affine  (Prandtl- 
Glauert)  space. 

The  design  methods  are  based  on  the  second-order  small-perturbation  theory 

4 

of  Weber  ; the  idea  is  to  obtain  an  approximation  for  the  planar  source  distri- 
bution, and  the  planar  doublet  distribution  when  appropriate,  compute  the  flow 
fields  due  to  these  singularities,  calculate  the  residual  errors  in  the  boundary 
conditions  on  the  upper  and  lower  surfaces,  and  adjust  the  source  distribution 
using  the  symmetrical  part  of  the  error  field,  and  the  camber  and  twist  distri- 
bution using  the  antisymmetrical  part.  Four  problems  are  studied. 

In  the  first  problem  the  thickness  and  doublet  distributions  are  prescribed; 

the  doublet  distribution  is  equivalent  to  the  load  distribution  in  linear  theory, 

and  is  roughly  equal  to  the  difference  in  lower  and  upper  surface  pressure 

coefficients  on  the  real,  thick  wing.  The  upper-surface  (and  of  course  the 

lower-surface)  pressure  distribution  is  found  as  part  of  the  solution,  along 

with  the  camber  and  twist.  This  method,  which  requires  nothing  beyond  the 

4 

suggestions  of  Weber  , converges  very  quickly,  and  takes  comparatively  little 
computer  time,  as  the  velocity  field  on  the  fixed  thickness  surface  due  to  the 
fixed  doublet  distribution  can  be  calculated  once  for  all. 

In  the  second  problem,  the  thickness  and  upper-surface  pressure  distribu- 
tions are  prescribed.  The  doublet  strength  now  has  to  be  determined  iteratively. 

4 

Weber  has  suggested  a procedure  for  this  task,  in  which  we  seek  to  cancel  the 
difference  between  the  target  pressure  distribution  and  that  achieved  so  far,  by 
a simple  addition  of  suitable  planar  doublet  strength  at  each  point,  the  idea 
being  that  the  resulting  extra  strearawash  will  dominate.  For  the  cases  tried, 
this  method  converges  quite  well  in  mid-semispan  and  - perhaps  surprisingly  - 
near  the  tip,  but  near  the  root  of  a swept  wing  there  are  difficulties,  in  that 
a small  change  in  doublet  strength,  particularly  near  the  root  trailing  edge, 
does  not  necessarily  produce  a small  change  in  the  other  two  components, 
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sidewash  and  upwash,  or  in  the  residual  error  which  determines  the  next  source 
distribution.  This  leads  to  a situation  of  the  tail  wagging  the  dog.  Various 
under-relaxation  schemes  have  been  tried,  with  limited  success.  More  successful 
was  an  inner  iteration  scheme  in  which  the  extra  sidewash  and  upwash  due  to  the 
new  doublet  increments  are  estimated  (on  the  chordal  surface,  with  a direct 
vortex  lattice  representation  for  the  upwash),  the  new  upper-surface  pressure 
distribution  estimated  and  another  increment  in  doublet  strength  obtained.  This 
inner  iteration  has  increased  the  scope  of  the  method,  but  there  are  still  some 
cases  when  it  fails  to  converge,  probably  because  for  the  prevailing  local  and 
temporal  conditions  the  pressure  is  bounded  away  from  the  target,  just  as  a 
quadratic  is  bounded  away  from  certain  regions.  For  this  case  an  optimization 
routine  would  be  appropriate,  but  has  not  been  included  in  the  program.  Thus 
at  the  moment,  this  option  cannot  be  guaranteed  to  work  near  the  root  of  a swept 
wing. 

In  the  third  problem,  the  loading  and  upper-surface  pressure  distributions 
are  prescribed,  and  the  thickness  is  to  be  determined  as  well  as  the  camber  and 
twist.  Since  it  is  the  doublet  strength,  and  not  the  exact  lower-surface 
pressure  distribution,  that  is  prescribed,  this  is  not  the  same  thing  as  the 
specification  of  both  upper-  and  lower-surface  pressure  distributions,  only 
nearly  so.  This  time,  after  the  computation  of  the  first  set  of  velocity  fields, 
the  shortfall  in  suction  is  regarded  as  a shortfall  in  streamwash  due  to  sources, 
and  this  is  turned  into  a perturbation  thickness  distribution  by  making  use  of 
the  approximate  formula  of  the  RAE  Standard  Method^  as  suggested  in  earlier  work 
by  Weber^.  Good  convergence  is  achieved  everywhere  except  near  the  root  of  a 
swept  wing,  where  over-relaxation,  or  an  inner  iteration  scheme  similar  to  that 
described  for  the  second  problem,  to  take  account  roughly  of  the  changes  in 
velocity  components  due  to  the  change  in  the  thickness  surface,  might  speed  up 
convergence  but  have  not  been  programmed. 

The  fourth  problem  is  a hybrid  of  the  second  and  third  problems.  Again 
the  upper-surface  pressure  distribution  is  prescribed  everywhere.  A chordwise 
section,  h = n*  in  the  non-dimensional  spanwise  variable,  is  chosen  near  the 
root.  Outboard  of  this  section,  the  thickness  distribution  is  prescribed,  as 
in  the  second  problem.  This  forms  the  basis  for  an  iterative  computation  of  the 
outboard  doublet  distribution.  The  condition  is  now  imposed  that  inboard  of 
q * q*  the  chordwise  variation  of  the  doublet  distribution  shall  be  the  same 

= n*  . In  this  way  we  expect  to  be  able  to  maintain  good  upper  surface 


as  at  n 
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flow  quality  and  wing  loading  right  into  the  root,  which  should  also  be  useful  in 
obtaining  root  stall  rather  than  tip  stall  at  high  off-design  incidence. 

Finally,  to  secure  the  required  upper-surface  pressure  distribution  inboard  of 
n = n*  , the  thickness  in  that  region  is  allowed  to  vary  as  in  the  third 
problem.  At  each  iteration  the  computed  thickness  distribution  is  faired  into 
the  given  outboard  distribution  by  a least  squares  method;  this  means  a certain 
loss  of  freedom,  as  the  upper-surface  pressure  condition  may  be  satisfied  only 
in  a spanwise  mean  sense,  but  it  is  judged  important  to  take  this  precaution,  in 
order  to  avoid  problems  in  final  manufacture  and  attachment  of  the  wing. 

As  for  the  second  problem,  good  convergence  is  found  in  mid-semispan  and 
near  the  tip,  and  we  can  expect  better  convergence  near  the  root  since  the  diffi- 
culty associated  with  the  second  problem  in  that  region  has  been  bypassed. 

Since  there  is  more  internal  freedom  than  in  the  third  problem,  the  most  we  can 
really  expect  a priori  is  that  convergence  near  the  root  will  not  be  worse  than 
for  that  problem,  and  in  fact  it  is  slightly  better  for  the  case  reported  here. 

In  our  example  cases  for  the  second  and  third  problems,  the  final  results 
are  noticeably  different  from  those  which  we  would  have  obtained  using  first- 
order  theory,  which  has  been  the  basis  of  classical  design  methods  for  many 
years,  and  is  used  here  to  provide  first  estimates  with  which  to  start  the 
iterations.  Thus,  these  programs  should  at  least  be  useful  tools  for  checking 
and  refining  the  results  of  classical  low-speed  wing  design  calculations.  It 
seems  likely  that  the  program  for  the  fourth  problem,  which  treats  what  seems  to 
be  a quite  practical  mixture  of  design  conditions,  will  also  be  useful. 

For  the  cases  considered  here,  which  are  slight  variations  of  RAE  wing 

8 

'B*  , acceptable  convergence  was  obtained  for  the  first  problem  in  two  or  three 
iterations,  for  the  second  in  four,  for  the  third  in  five  (except  near  the  root) 
and  for  the  fourth  in  five  (including  the  root).  Thus  the  second,  third  and 
fourth  problems  take  around  twice  as  long  as  the  direct  program*.  This  is 

unfortunate,  and  may  mean  that  the  method  is  less  competitive  than  the  repeated 

9 ... 

use  of  the  BAC  program  in  an  incremental  mode.  But  there  is  the  possibility  of 

performing  only  one  iteration  at  small  cost  and  having  a look  at  the  results, 

then  deciding  whether  to  continue  or  to  change  an  unpromising  set  of  design 

distributions  and  start  again,  without  incurring  a large  penalty  in  computer 

time. 
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I Two  separate  programs  have  been  written,  one  for  the  first  and  second 

problems,  one  for  the  third  and  fourth.  The  programs  are  written  in  Fortran 
and  occupy  about  45K  words  of  core  store,  with  arrays  dimensioned  for  a 12  x 15 
collocation  grid.  The  vortex  lattice  matrix  is  stored  on  a disc. 

12  BASIC  EQUATIONS 

In  this  section,  we  summarize  the  basic  equations  and  boundary  conditions 
as  set  out  in  detail  in  Ref.l. 

We  consider  a finite  wing  in  a uniform  free  stream  with  speed  unity  and 

Mach  number  < 1 . We  take  cartesian  coordinates  (x*,y,z*)  with  origin  0* 

at  the  apex  of  the  wing,  with  the  x*-axis  0*x*  in  the  free  stream  direction, 

0*y  to  starboard  and  0*z*  upwards,  z*  = 0 then  defining  the  datum  plane  of 

the  wing,  so  that  at  any  plane  station  y = constant,  the  local  wing  section 

incidence  is  just  the  angle  of  twist  a^,(y)  ; in  what  follows  the  distribution 

has  to  be  determined  as  part  of  the  design  problem.  In  any  plane 

y = constant,  we  define  local  cartesian  coordinates  (x,z)  such  that  the  x-axis 

is  parallel  to  the  local  section  chordline,  the  z-axis  completes  a right-handed 

set  with  the  x-  and  y-axes,  and  the  origin  0 is  the  z-projection  of  the  apex  0* 

(Fig.l).  These  axes  can  be  got  by  rotating  the  (x*,y*)-axes  through  the  angle 

a^,(y)  . Wing  thickness  (z^)  and  camber  (zg)  ordinates  are  defined  as  ordinates 

z normal  to  the  local  chord: 
w 

z„(x,y)  = ± z (x,y)  + z (x,y)  . (2-1) 

W L s 

Both  zt  and  zg  vanish  at  leading  and  trailing  edges. 

The  velocity  in  the  (x,y,z)  system  is  now  made  up  of  the  free  stream 
velocity 

U = (cos  ctT,0,sin  aT>  (2-2) 

and  a perturbation  velocity  field 

u = (u,v,w) 

which  is  to  satisfy  the  wing  surface  boundary  condition 

(U^  + u)  . grad(z^  - z)  =0  . (2-3) 

For  compressible  flow,  with  > 0 , we  transform  into  the  affine 
(Prandtl-Glauert , or  Goethert)  space  (x,y,z)  with  the  corresponding  affine 
perturbation  velocity 
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r 


u = (u,v,w) 


by 


(2-4) 


where  6“ 


1 - M 


To  a first  approximation,  the  problem  is  reduced  to  an  incompressible  flow 

problem  in  the  affine  space,  which  we  proceed  to  tackle  by  placing  source  and 

doublet  distributions  q(x,y),  £(x,y)  on  the  chordal  surface  of  the  analogous 

wing.  These  distributions  will  generate  the  affine  perturbation  velocity  field 

u . We  assume  that  at  each  station  y = constant  , u will  not  be  significantly 

affected  if  the  source  and  doublet  distributions  are  considered  as  lying  in  the 

2 3 

local  plane  z = 0 , so  that  it  can  be  calculated  using  the  Ledger-Sells  ’ 
computer  subroutine. 


Using  this  subroutine,  the  separate  perturbation  velocities  ii^,  Ci  , due 
to  sources  and  doublets  respectively,  will  be  calculated  on  the  'thickness 
surface ' 

derived  using  the  first  two  terms  of  a Taylor  expansion  in  z 


z = z , and  their  values  on  the  actual  wing  surface  z = z will  be 
t w 

The  upwash 

0 

w£  due  to  the  sources,  and  the  streamwash  u^  and  sidewash  v^  due  to  the 
doublets,  change  sign  on  going  from  z = z^  to  z = -z  , and  so  the  components 
are : 


of  u 
— t 


and  u 
—l 


u 

— t 


L + 3Gt 

u ± z 

t s 9z 


Sv, 

v ± z 

t s Jz 


3w 

± w + z - — - 

t s 8z 


( , 3Q« 

r + zs  3z 


3v« 

± v2  + Zs  3z 


3w 

S 

w + z - — 
£ s 3 z 


evaluated  at  z = . Upper  and  lower  signs  correspond  to  upper  and  lower 

wing  surfaces,  respectively. 

For  the  complete  velocity  U=U  +u  +n=  (U,V,W)  , making  use  of 
(2-4),  we  then  have 


i / 8M  / 8at 

cos  a + - [ U + Z ± Z -r— - + U 

T S It  s 32  s 3z  £, 


v + z 
t s 3z 


s 3z  A 


(2-5) 


/ 3wt  \ / 8w 

W = sin  a + [z  — £ + w„  ± w + z — - 
T s 3z  l]  t s 3z 


These  equations  can  now  be  used  in  the  boundary  condition  (2-3)  and  will 

also  be  needed  to  obtain  the  total  velocity  Q from  Q2  = U2  + V2  + W2  and  the 

pressure  coefficient  C from 

P 


C = 1 - O'1 

P 


(M  = 0) 


[l  + J(Y  - !)M2(!  - Q2)] 


Y/(Y-1) 


(M  > 0)  , 

oo  * 


where  Y is  the  ratio  of  specific  heats,  taken  here  as  1.4. 

Within  second-order  theory,  for  the  boundary  condition  and  for  the  calcula- 
tion of  0^  we  can  write  sin  a^,  = in  (2-7),  and  for  the  boundary  condition 
(but  not  when  calculating  C^)  we  can  write  cos  «T  = 1 in  (2-5).  Making  these 
changes , we  have 


u = Q0  ± Q, 
v = q2±q3 
w = q5  ± q4 


r 


A 


(~  3Gi 

1 + ut  + zs  IT 


sir  ♦ 0 A 


♦ ,v« 
V + 2 

t s 3z 


3v 

z ^ — + v„ 

s 5z  2. 


3w 

w + z — ^ 
t s 3z 


3w 

t ^ 

a + z ■-  ■ + w 
T s 3 z l 


> 


(2-10) 


Substituting  from  (2-1)  and  (2-9)  in  the  boundary  condition  (2-3),  we  find 


± Rt  + Re  = 0 


whence  we  obtain  the  symmetric  part  of  the  boundary  condition 


R = 0 

t 


(2-11) 


and  the  antisymmetric  part 


R?,  = o 


residuals 

R^  and 

h 

are  given 

by 

R 

t 

%d_\ 

9z 

■IT  *q2 

32t 

+ Q3 

3z 

s 

6 3* 

6 

3y 

3y 

R* 

% 3Zs 

+ ^1 

S2t  * „ 

3z 

s 

3zt 

= T3T 

8 

' 3x  + Q2 

3y 

+ Q3 

3y 

- Qy 


- Q 


5 ’ 


(2-12) 


(2-13) 


(2-14) 
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It  is  convenient  here  to  introduce  the  coordinate  system  used  for  compu- 
tations, which  is  the  same  as  that  of  Ref.l.  The  local  percentage  chord  is 
given  by 

x = xL(y)  + c(y)£  (2-15) 

where  x^(y)  is  the  leading-edge  ordinate,  and  c(y)  the  local  chord,  of  the 
analogous  wing.  The  computation  grid  is  then  defined  chordwise  by  equal 
intervals  of  the  angular  chordwise  coordinate  $ where 

5 = iO  - cos  <f>)  . 

We  use  the  nondimens ional  spanwise  coordinate  n = y/s  . (The  semispan  s will 
be  taken  as  1 . ) We  also  define  a spanwise  variable  q = y to  go  with  £ , and 
use  d/dri  to  denote  partial  derivative  with  respect  to  fi  along  lines  of 
constant  percentage  chord  £ . Then 


JL  - JL  - an  A — 

3y  37]  9x 


/ L dc  \ 

where  A = arctan  I + — Cj  is  the  local  sweep  angle  (for  the  analogous  wing), 
3 THE  DESIGN  PROBLEMS  STUDIED 

3. 1 First  problem:  specified  loading  and  thickness  distributions 


In  this  problem,  the  camber  and  twist  distributions  zg,  are  to  be 

determined,  and  the  upper- surface  pressure  distribution  is  found  as  part  of  the 
solution.  We  have  only  to  determine,  in  addition,  the  auxiliary  unknown  source 
distribution,  since  the  doublet  distribution  J,(x,y)  (which  is  equivalent  to 
the  loading  distribution  in  linear  theory)  is  known  in  advance. 

Since  the  doublet  distribution  is  fixed,  we  can  compute  the  velocity 

field  u on  the  thickness  surface  z = z once  for  all.  Let  us  assume  that 

“«•  • c(n) 
we  have  also  determined  an  approximation  q (x,y)  to  the  final  source 

distribution,  and  that  we  have  computed  the  corresponding  velocity  field  . 


and  twist 


We  suppose  also  that  estimates  for  the  camber  distribution  z^ 

/ \ s 

distribution  a,:  are  available.  Using  (2-13),  (2-14),  we  compute  the 

l / \ / \ 

corresponding  residuals  R w ' and  examine  how  well  the  boundary  conditions 


(2-11),  (2-12)  are  satisfied.  Consider  the  first  of  these.  The  residual  R^ 
can  be  thought  of  as  a deficiency  in  w through  the  tern  (-Q,),  and  so  we 

, *4 


(n) 


attempt  to  cancel  it  by  adding  to  q^n^ 

Aq  = 2R 


a source  distribution 
(n) 


(3-1) 


exactly  as  in  Ref.l.  To  start  the  iteration  and  obtain  the  first  estimate 
(the  basic  source  distribution  q^) , we  take  the  situation  where  0^^  , ^ ^ 
are  both  zero,  giving  the  linear-theory  result 


,(0) 


llll 

e 3x 


3z 

t 

3x 


again  as  in  Ref. 1 . 

The  residual  error  R^ 
adjust  the  current  values  z^n^ , a,£n)  of  z^ , as  shown  by  WeberA.  Putting 


(n) 


in  the  other  boundary  condition  can  be  used  to 


z = z(n)  + Az 
s s s 


4n>  * ‘“i 


in  (2-14),  neglecting  products  of  the  perturbations  Az  , Aa^,  with  other 
perturbation  quantities,  and  invoking  (2-12),  we  have 


. 3 Az 

1 s 

7 3x  AaT 


= - R 


(n) 


(3-2) 


ZsB*  aTB 


with  which  to  start 


To  obtain  the  basic  camber  and  twist  distributions 
the  iteration,  we  compute  W£B^>y*zt^  or  esti1113116  from  the  basic 

doublet  distribution 


l , and  then  (2-14)  gives 

D 


,(0) 


£B 


to  be  substituted  in  (3-2),  which  leads  to  the  standard  result  of  linear  wing 
design  theory. 


The  integrals  in  (3-3),  (3-4)  are  evaluated  using  Simpson's  rule  with  <(> 

as  the  independent  variable.  Then  the  accumulated  z . a„,  are  extrapolated  to 

s ’ T 

the  wing  root  and  tip  stations,  and  the  spanwise  derivative  9zg/3ri  (also 
needed  in  the  Taylor  series  expansions)  is  computed  at  each  collocation  point 
by  a cubic  spline  fitting  routine. 

We  remark  that  the  calculation  of  Az  , Act_  does  not  disturb  the  field 

s T 

values  u^,  on  z = , and  so  after  finding  Azg,  AoT  via  the  residual 

error  field  from  Taylor  series,  it  is  worthwhile  to  execute  the  Taylor  series 
sequence  again  with  the  updated  values  of  z , a . This  updates  the  successive 

S i 

estimates  for  the  other  residual  field  , as  well  as  the  upper-surface 

pressure  distribution,  which  we  shall  need  in  the  other  design  problems  to  be 
considered. 

3.2  Second  problem:  specified  thickness  and  upper-surface  pressure 
distributions 


In  this  problem,  the  doublet  distribution  H is  to  be  determined 
iteratively,  as  well  as  the  source  distribution  q and  the  camber  and  twist 
distributions  z , a . A sequence  of  upper-surface  pressure  distributions 

is  generated,  which  (it  is  hoped)  will  tend  to  the  specified  or  target 

^ ...  - . . 4 . . 

distribution  C . Again  following  Weber  in  broad  outline,  successive 

increments  AH  in  H are  determined  from  the  shortfall  in  C ; then,  iust 

pu  ’ 

as  in  the  first  problem,  successive  increments  Aq  arc  found  from  the  sequence 

of  residual  errors  , and  successive  increments  Az  , Aa_  from  the 

t ( n ) s X 

sequence  of  residual  errors  R^ 

A first  estimate  H for  i.  is  provided  by  Lock^.  From  the  target 

D 

distribution  C and  the  appropriate  equation  (2-8),  we  derive  the  target 
upper-surface  total  speed  distribution  Q , and  insert  it  into  Lock's  formula: 


I 
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DQ2 


[COS  “t  + UtB  C°S  °T  + Ul  B ('  + S(3>  S6C  AS/6n)j 

+ [vtB  C°S  °T  + v«(*  + s0)  SCC  Am/6n)] 

/ 2\  2 
+ ^1  - k2J(D  - s^n 


‘x  cos2“t 


(3-5) 


The  symbols  D , etc.  are  defined  in  Appendix  A.  The  velocity  components  are 
taken  in  the  plane  z = 0 , so  that  the  doublet  strength  is  connected  to  the 
(physical)  component  u _ by 

X.O 


eB  A6uiB  - 


If  we  tentatively  assume  that  v ^ 
relation 


depends  linearly  on 


JIB 


through  a 


v 


£B 


U*B  Can  Av 


(3-6) 


where  is  also  defined  in  Appendix  A,  then  equation  (3-5)  becomes  a quad- 

ratic in  u^B  (all  other  quantities  being  known),  and  one  of  the  two  roots  can 
be  picked  out  as  the  required  solution. 

We  can  improve  this  first  estimate  in  an  iterative  way,  calculating  v^ 
from  the  current  estimate  for  u,  and  the  equation  of  irrotational  flow  and 
using  this  value  instead  of  the  assumption  (3-6).  Details  are  again  set  out  in 
Appendix  A.  This  is  in  principle  an  inner  iteration,  performed  before  the 
velocity  fields  u£*\  u^* ^ are  computed  accurately  in  the  main  iteration 
cycle.  In  practice,  one  inner  iteration  is  sufficient. 

In  an  inner  iteration  cycle  such  as  this,  in  order  to  estimate  the 
successive  camber  and  twist  distributions  ve  would  like  to  obtain  a quick 
estimate  for  the  upwash  corresponding  to  a current,  or  intermediate,  estimate  of 
the  doublet  field  without  performing  the  time-consuming  Ledger-Sells  double 
integration  for  each  intermediate  estimate  in  turn.  To  do  this,  we  use  a direct 
vortex  lattice  representation,  as  in  Ref.l;  the  vortex  lattice  influence  matrix 
is  not  inverted  (but  as  appreciable  time  is  needed  to  generate  it,  it  is  stored 
on  a scratch  disc  by  the  computer  after  the  first  pass  through  the  vortex 
lattice  subroutine).  Near  the  root  and  tip,  the  vortex  lattice  method  is  not 
very  accurate  and  tends  to  overestimate  the  upwash  o’  .ownwash,  and  so  the  out- 
put values  are  multiplied  by  a spanwise  under-relaxation  factor  tentatively 

taken  (after  numerical  experiments)  as: 
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db  - 0 

V 0 + l<2l) 

k-,  is  a spanwise  interpolation  function,  similar  to  that  of  the  RAE  Standard 
Method,  and  defined  in  Ref.l.  The  errors  in  the  resulting  estimates  for  camber 
and  twist  will  not  matter  greatly,  since  they  will  be  nearly  corrected  auto- 
matically when  u^  is  computed  accurately  in  the  main  iteration  cycle. 

As  remarked  earlier,  after  the  calculation  of  and  at  the 

Pu 

end  of  the  nth  main  iteration,  we  have  next  to  generate  an  increment  AH  in  the 

doublet  strength  from  the  shortfall  in  C . This  shortfall  is,  in  principle, 

Pu  4 

a second-order  quantity  and  from  it,  following  Weber  , we  can  derive  a second- 
order  expression  for  the  corresponding  increment  Au^  in  (physical)  streamwash 
due  to  AH  , where  as  usual  AH  = 4f3Auf  . Using  the  suffix  u to  represent 
upper-surface  values  in  the  velocity  components  (U,V,W)  from  (2-5)  to  (2-7),  we 

, v 2 2 2 2 

have  Q'  =U  + V + W , and  to  second-order  accuracy 
u u u u 

_2  2 2 2 
Q = (U  + Au. ) + V + W 

u H u u 

= 0(n)2  + 2U  Au 

'u  u H 


whence 


Q2  - Q(n)2 

u 


We  now  have  all  the  equations  needed  to  set  up  a closed  iteration  cycle  for 
this  problem. 

In  order  to  control  some  overshoot  near  the  tip  (and  as  a partial  control 
near  the  root  also  - this  will  be  discussed  further,  later),  it  has  been  found 
helpful  to  introduce  an  under-relaxation  factor  for  the  doublet  strength 

perturbations,  again  depending  on  the  spanwise  factor  <2  : 


% = o * l<2|)  5 . 

Near  a wing  tip  of  finite  chord,  the  loading  is  expected  (apart  from  the  effect 
of  corner  singularities)  to  decay  elliptically ; the  factor  is  not  expected 

to  represent  this  spanwise  decay  precisely,  but  is  intended  to  introduce  some 
decay,  under  control,  which  seems  in  practice  to  be  better  than  introducing  none 
at  all. 


15 


It  will  be  remembered*  that,  to  assist  the  solution  of  the  direct  problem, 
line  source  and  doublet  distributions  were  introduced  on  the  wing  centre  line 
y = z = 0 to  reduce  the  residual  errors  Rt,  Rf  at  the  root  y = 0 . In  the 
present  design  problem,  just  as  we  adjust  the  planar  doublet  strength  to  produce 
an  extra  Au^  to  correct  the  upper-surface  pressure  distribution  outboard,  it 
is  natural  to  try  to  adjust  the  line  doublet  strength  to  produce  an  extra  stream- 
wash  Au^  to  do  the  same  job  at  the  root.  The  resulting  upwash  adjustment 
AW£D  would  then  be  incorporated  into  the  camber  and  twist  distribution  at  and 
near  the  root.  The  difficulty  is,  that  whereas  (in  the  direct  problem)  a change 
in  line  doublet  strength  to  produce  a certain  Awf^  did  not  tend  to  produce  a 
large  change  Au^  » a change  to  produce  a certain  Au^  does  tend  to  produce  a 
large  change  Aw^  ; in  other  words,  this  calculation,  although  well-conditioned 
in  one  direction,  is  ill-conditioned  in  the  other.  It  was  found  that  line 
sources  also  tended  to  destabilize  the  iteration  scheme;  so  the  line  singularity 
technique  has  been  abandoned  and  no  target  distribution  is  prescribed  at  the 
wing  root  (though  we  can  approach  it  quite  closely  with  a suitable  choice  of 
spanwise  stations);  all  quantities  such  as  z , q , i are  extrapolated  para- 
bolically  to  the  root  (but  , u^  are  still  evaluated  using  the  Ledger-Sells 
subroutine),  and  the  final  residual  errors  and  pressure  coefficients  at  the  root 
are  left  to  take  care  of  themselves. 

As  in  the  direct  program* , we  can  seek  to  reduce  the  number  of  iterations 
by  generating  improved  estimates  for  the  perturbation  quantities.  The  additional 
velocity  fields  due  to  the  perturbation  source  and  doublet  fields,  calculated 
in  the  main  iteration  cycle,  can  be  estimated  using  the  RAE  Standard  Method^  and 
so  it  is  possible  in  effect  to  perform  one  iteration  cycle  and  to  obtain  a 
further  set  of  perturbation  quantities,  without  actually  performing  the  accurate 
but  lengthy  Ledger-Sells  calculations.  To  estimate  the  new  residual  errors,  we 
again  invoke  Maclaurin  series  (expansions  about  z = 0 , rather  than  z = zt)» 
modified  to  take  account  of  concurrent  changes  in  camber.  The  algebra  is  set 
out  in  Appendix  B.  This  can  indeed  be  incorporated  in  a further  inner  iteration 
cycle,  but  it  has  been  found  inadvisable  to  do  it  more  than  once  or  twice; 
convergence  of  the  main  iteration  cycle  is  rather  sensitive,  perhaps  because  of 
the  level  of  feedback  involved  between  the  three  unknown  and  interacting  field 
quantities,  and  it  seems  likely  that  any  further  benefit  derived  after  the  first 
or  second  inner  iteration  is  outweighed  by  the  feedback  effect  due  to  th 
accumulating  errors  in  the  successive  estimates  from  each  cycle. 


- 
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In  the  course  of  developing  the  program,  a difficulty  has  arisen  which  has 

still  not  been  satisfactorily  resolved.  The  symptom  of  this  difficulty  is  the 

onset  of  oscillations  or  divergence  in  v„ , w . C and  R near  the  root 

° Z Z pu  t 

trailing  edge.  In  this  region  there  seems  to  be  a close  coupling  between  all 
the  velocity  components  so  that  a small  change  in  l and  hence  in  u^  does  not 


necessarily  produce  a negligible  change  in  v^,  w^ 


; thus  the  effect  of  the  small  change 
required  effect  on  Q 


AL 


and  in  the  residual  error 
given  by  (3-7)  does  not  have  the 


To  prevent  R and  the  resulting  extra  source  distributions  from  growing 
too  large,  after  the  first  guesses  for  source  and  doublet  fields  the  program  has 
been  arranged  to  do  two  successive  perturbation  source  calculations  (including 
Maclaurin  series  calculations),  just  as  if  it  were  performing  two  iterations  of 
the  first  problem,  before  calculating  another  perturbation  doublet  field.  This 
is  an  unfortunate  necessity,  but  as  the  Ledger-Sells  subroutine  is  adaptive,  the 
calculation  times  should  decrease  as  the  iterations  proceed  and  R^  decreases. 


The  destabilizing  effect  of  Av^  and  Aw^  is  not  so  easily  dealt  with. 

An  attempt  to  take  them  into  account  in  the  derivation  of  equation  (3-7),  even 
whep  only  linear  terms  are  included,  leads  to  some  heavy  matrix  algebra  and 
programing,  and  fails  in  the  end  because  the  non-linear  terms  are  not  negligible 
in  practice.  For  some  of  the  cases  studied,  instability  was  successfully 
averted  by  a simple  addition  to  the  inner  iteration  scheme,  in  which  we  take 
account  of  the  approximate  values  of  AQ^,  Av^,  Aw^  (and  also  of  Aut,  Avt , Aw^ 
generated  within  the  Maclaurin  series  sequence,  when  appropriate)  to  update 
approximately  the  upper-surface  velocity  components  U^,  V , and  the  total 

speed  and  hence  to  generate  a further  perturbation  doublet  field  using  (3-7) 

again.  This  part  of  the  inner  iteration  cycle  can  now  be  repeated  until  the 
changes  in  doublet  strength  are  sufficiently  small,  or  for  a maximum  of 
(currently)  nine  inner  cycles.  This  inner  iteration  scheme  can  also  be 
profitably  applied  to  the  rest  of  the  wing  surface,  even  though  the  effect  of 

Av„  and  Aw„  on  Q is  not  so  great. 

Z Z u 6 

This  artifice  has  enabled  us  to  obtain  a solution  for  at  least  one  case 
which  we  could  not  treat  without  it.  But  there  seems  to  be  a class  of  cases  for 
which  the  program  so  far  described  still  does  not  work;  for  these  cases,  it  is 
the  inner  iteration  which  oscillates  or  diverges.  Relaxation  methods  do  not 
seem  to  help.  A study  of  the  computer  output  indicates  that,  even  though  a 
numerically  exact  solution  of  the  whole  problem  is  known  (derived  from  the 
program  for  the  first  problem,  for  instance),  the  unconverged  state  of  the 
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accurately  calculated  parts  of  u^  and  IS  , and  the  approximations  used  for 

AG.  in  terms  of  M in  the  inner  iteration  scheme,  are  such  that  Q has  a 
— r u 

minimum  value  which  is  higher  than  the  design  value  Q , just  as  a quadratic 
in  a real  variable  is  bounded  away  from  certain  values.  Thus  in  this  situation 
we  cannot  find  a M to  give  Qy  = Q everywhere,  using  just  the  inner 
iteration  scheme.  An  optimization  sequence  within  this  scheme,  to  detect  the 
situation  and  to  seek  a solution  minimizing  the  overall  differences  between  Q 
and  Q (by  least  squares,  perhaps)  is  indicated,  but  has  not  been  devised  at 
the  time  of  writing. 


3. 3 Third  problem:  specified  loading  and  upper-surface  pressure  distributions 

In  this  problem,  the  thickness  distribution  z^_  is  to  be  determined 
iteratively,  as  well  as  the  camber  and  twist  distributions  z 
associated  source  distribution  q 


, (n) 


s'  T 

Again  a sequence  of  upper-surface  pressure 


distributions  C^“/  is  generated,  which  (it  is  hoped)  will  tend  to  the  target 
distribution  C 

pu 

As  in  the  first  and  second  problems,  successive  source  increments  Aq  are 

(n) 


found  from  the  sequence  of  residual  errors 
twist  increments  Az 


and  successive  camber  and 
(n) 


Act^  from  the  sequence  of  residual  errors 


R. 


This 


leaves  the  successive  thickness  increments 


Az^  to  be  determined  from  the 


shortfall  in  upper-surface  C 


pu 


We  can  do  this  approximately,  making  use  of 


results  from  the  RAE  Standard  Method  , if  this  shortfall  can  be  converted 


approximately  into  a shortfall 
three  sets  of  circumstances: 


Au. 


in  u 


We  may  consider  the  following 


(i)  At  the  outset,  we  have  no  estimate  at  all  for  z^  . Given  the 
design  upper-surface  velocity  and  the  chordal-surface  streamwash  u^g  due 

to  the  specified  load  (doublet)  distribution,  we  have  a simple  basic  estimate 
for  the  streamwash  due  to  sources: 


atB  = *(QU  " 0 " S*B  * 

This  will  lead  to  a first  estimate  z „ for  the  thickness  distribution. 

tB 

(ii)  We  can  modify  this  initial  crude  estimate  for  z with  the  help  of 
Lock's  formula  (3-3).  From  the  known  doublet  strength,  we  again  have  upR 
(hence  u .„)  and  also  v.„  , on  the  chordal  surface.  If  we  assume  that  u._  , 

)CD  1 CD  tr 

in  the  first  term  on  the  right  side  of  (3-5),  dominates  the  equation  as  far  as 
sources  are  concerned,  we  can  substitute  for  all  the  other  source  terms,  such  as 


18 


v _ , the  values  determined  from  the  initial  estimate  z _ and  obtain  a 
tB  * tB 

revised  estimate  for  u (and  hence  a revised  estimate  for  z ) . The  details 

t t 

are  filled  in  at  the  end  of  Appendix  A. 

(iii)  After  computing  in  the  nth  main  iteration,  just  as  in  the 

second  problem  we  can  derive  a simple  second-order  expression  for  the  required 
physical  increment  Aut  = Au^/8  due  to  the  required  extra  thickness  distribution 
Az  . This  expression  is  just  the  right  side  of  equation  (3-7)  again: 

4“t  ■ - Qin)2]/2U u • 

After  a Maclaurin  series  cycle,  the  velocity  components  in  this  expression  can 
be  approximately  updated,  as  in  the  second  problem. 

The  perturbation  velocity  Afl^  and  the  extra  thickness  Az^  (which  can 

be  thought  of  as  an  extra  source  distribution  23Az  /3x)  are  connected  approxi- 

£ L 

mately  by  the  formula  of  the  RAE  Standard  Method  : 


/9Az  \ _ 3Az  /3x 

Au  = cos  A a(  ) - x0(n)f(A)  ; ^ 

C \ / Jl  ♦ (3Azt/3*)‘ 


1 f /9Azt 
0 7T  J ( 3S 


dt' 

~ T 


. 1,1+  sin  A 

f = — In  — - 

11  1 + sin  A 


<2  is  the  spanwise  interpolation  function,  defined  in  Ref.l.  One  way  of 
dealing  with  this  equation  is  to  set  up  an  iterative  cycle  in  which  the  last 
term  (with  ’Riegels  factor')  is  treated  as  small  and  known  from  the  previous 
iteration  - indeed,  it  vanishes  altogether  at  certain  mid-wing  stations  where 


<2=0.  Thus  we  write 


/3Az(m) 

t 


+ <-f(A) 

cos  A z 


fl  + (3Az 


5Az^m_ ' ^ I 3xj‘ 


(3-9) 


at  the  mth  such  cycle,  and  to  start  the  cycle  we  put  Az^^  = 0 . The  solution 
of  this  equation  (see  Ref. 7,  for  example)  with  Aztm  (1)  = Az^  (0)  = 0 is: 


T 
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9Az 


(m) 


3x 


V c (i  - 6 


7 / - <’>  r^V 


(3-10) 


The  program  already  incorporates  a subroutine  to  evaluate  the  Cauchy 
(principal-value)  integral  (3-8),  and  the  Cauchy  integral  (3-10)  can  be  evaluated 
using  the  same  subroutine;  thus  the  method  is  more  convenient  than,  for  example, 
an  extension  of  Carleman’s  method  as  in  Ref. 7. 

However,  before  actually  coding  the  evaluation  of  (3-10)  we  should  consider 

that  the  initial  guess  z _ , and  the  corrections  Az  , should  represent  a 

tb  t 

wedge-shaped  trailing-edge  (£  = 1)  for  structural  reasons;  whereas,  if  F is 
assumed  regular  at  E,  = 1 , then  9Azt/9x  is  o[(i  - and  Az^_  takes  a 

rounded  or  elliptical  shape  there,  as  well  as  at  the  leading-edge  £ = 0 . On 
the  other  hand,  when  Az^_  is  regular,  a and  F will  show  a logarithmic 
singularity  as  £ -*•  1 . We  therefore  estimate  the  strength  of  the  singularity 
and  subtract  a suitable  function  from  F to  leave  an  integrand  which  we  hope 
will  lead  to  a sensible  trailing-edge  shape. 

A representative  function  A^  which  represents  a closed  section,  with  a 
suitably  small  square-root  singularity  at  the  leading-edge  and  wedge-shaped  at 
the  trailing-edge,  is 


3z.-n 

1 = AQ  = 4(3  — 5£)/£  (3-)l) 

which  corresponds  to  a section  shape 

zt0/A  = (I  - OK2'2  . 


Inserting  (3-11)  in  (3-8),  we  have 


a 


0 


2 

7T 


1 ♦ fj 
1 - /£ 


1 


We  now  subtract  a suitable  multiple  6 of  from  F . To  help 

determine  this  multiple,  we  have  the  values  of  F at  discrete  chordwise 
stations  E,^,  £^,  ...,  £^_  ^ , £^  where  0 = £ < ^ < • • • < ' • We  make 

the  hypothesis 
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F = SoQ(0  + Oj(5) 


where  0j(O  varies  slowly  compared  with  °q(0  near 
£ = CL_j,  in  turn  and  subtracting,  we  find 


£ = 1 . Putting 


6 = 


r(EL) 


We  then  evaluate  A*  from  the  equation,  similar  to  (3-10) 


A*  = 


VC(1  “ O 


] 

7 / [f(C')  - 6o0U')j^’(i  - S’) 


dC 


and  then 


3Az 


(m) 


355 


= A*  + 5A, 


This  can  now  be  substituted  into  F given  by  (3-9),  with  m increased  by  one, 
to  start  the  next  iteration  cycle.  We  have  tested  this  cycle  on  a variety  of 
wings  and  shapes,  and  convergence  always  seems  rapid,  up  to  10  iterations  being 
needed  depending  on  the  value  of  the  interpolation  function  < ^ . 

Although  the  integral  (3-10)  and  the  function  (3-11)  always  represent 
closed  contours  mathematically,  when  Az^ra^  is  evaluated  by  numerical  integra- 
tion of  3Az^m^ 351  this  may  not  be  exactly  true  because  of  numerical  truncation 
errors.  At  each  cycle,  then,  we  check  the  closure  condition  and  rotate  the 

contour  Az ^ through  a small  correction  angle  about  the  leading-edge  in  the 
C / \ 

(5,Azt)  plane  to  make  Az£  ' = 0 at  £ = 1 . 

After  calculating  Az^_  we  must  check  that  the  resulting  contour  does  not 
cross  itself,  i.e.  the  new  z^  is  positive  everywhere.  Actually,  we  demand 
rather  more.  To  avoid  the  possibility  of  unrealistically  small  z^  , we  check 
that  z^_  will  not  be  reduced  by  more  than  half  its  former  value  anywhere;  if 
this  condition  is  not  met,  i.e.  Azt  < - Jz^  , then  Az^  is  everywhere  multi- 
plied by  a suitable  factor  to  give  just  the  50  per  cent  reduction  in  z ^ at 
some  point  which  is  the  largest  reduction  we  are  prepared  to  allow. 
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On  the  other  hand,  if  the  first  estimate  z fails  this  test,  since  there 

is  no  previous  distribution  to  compare  it  with,  we  have  to  improvise.  First,  we 

examine  z near  the  leading  edge;  if  it  is  negative  somewhere  in  that  region, 

then  probably  the  input  data  for  l and  C do  not  correspond  to  a reasonable 

wing;  the  program  exits.  If  z is  positive  there,  but  becomes  negative  near 

the  wing  trailing-edge , we  charitably  assume  that  the  input  data  do  correspond 

to  a reasonable  wing  and  the  initial  guess  for  z is  in  error.  At  each  span- 

tB 

wise  station  we  then  find  the  maximum  z .at  E = E say.  and  also  note  the 

tB  * 

derivative  9z  / 3 E at  E . For  E < E < 1 we  construct  a function  G(E) 
tB 

which  is  arcwise  continuous  with  z at  E = E * becomes  rapidly  linear  as  we 

t B 

recede  from  that  point,  and  vanishes  at  the  trailing-edge  E = 1 . A suitable 
function  is 

G = ae"k(S_^  + bE  + c 
where  we  arbitrarily  choose  k so  that 

k(l  -0-2 


so  that  the  exponential  term  dies  out  near  E = 1 • a,  b and  c can  now  be 
chosen  to  meet  the  three  conditions  (Appendix  C). 

We  apply  this  procedure  at  all  spanwise  stations,  even  if  a negative 
ztg  has  been  detected  at  only  one  station,  because  9zt^/9y  (which  is  needed 
for  the  boundary  conditions)  has  to  be  computed  from  the  values  generated,  and 
if  we  only  adjusted  the  values  at  one  such  station,  the  difference  between  the 
adjusted  section  and  the  neighbouring  unadjusted  sections  downstream  of  E = E 
might  cause  considerable  fluctuations  in  this  spanwise  derivative.  By  applying 
the  procedure  uniformly  at  all  spanwise  stations,  we  hope  to  avoid  this  possible 
source  of  trouble. 

Indeed,  each  time  a fresh  z^  is  computed  for  this  problem,  we  have  to 

compute  and  store  the  new  chordwise  derivative  9z  / 9 E as  well  as  9z  /9y  , and 

1 *■ 
we  also  need  the  tables  of  arclengths  along  chordwise  and  spanwise  curves  on 

the  new  thickness  surface  z = z , and  the  Lighthill  E-shift  factors  rendering 

the  solution  uniformly  valid  near  the  new  rounded  leading-edge. 

Further,  the  device  of  computing  only  the  (small)  perturbation  velocities 
Au^,  Au^  from  Ae^,  Ae^  , which  saves  time  because  the  Ledger-Sells  routine  is 
adaptive  (useful  for  the  first  and  second  problems),  breaks  down  here  because 
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the  already-computed  velocities  are  known  only  on  some  previously-generated 
thickness  surface,  not  at  the  current  values  of  z . (An  attempt  to  calculate 
values  on  the  new  thickness  surface  using  Taylor  series  turned  out,  not  surpris- 
ingly, very  inaccurate  near  the  leading-edge.)  So  for  this  problem  it  is 
necessary  to  work  with  the  complete  source  and  doublet  distributions  all  the 
time  when  calculating  u^  and  u^  . 

When  the  program  so  far  described  was  run,  inspection  of  the  results 
showed  that  despite  the  inclusion  of  Az^  in  the  Maclaurin  series  theory,  the 
successive  estimates  for  Az^  would  be  improved  if  Az^_  was  determined  only 
after  one  Maclaurin  series  calculation,  and  not  after  the  Taylor  series. 

Similarly  to  the  second  problem,  there  seems  to  be  a coupling  between  R^,  Az^ 
and  the  shortfall  in  , such  that  if  Az^_  is  calculated  every  time,  even 

though  the  corresponding  linear-theory  source  distribution  is  always  taken 
into  account,  R^  does  not  converge  quickly  to  zero.  By  determining  Az^_ 
only  half  as  often,  the  residual  R^_  is  given  more  time  to  settle  down;  and  by 
determining  it  after  the  Maclaurin  series  calculations,  we  postpone  it  as  long 
as  possible  after  the  final  Ledger-Sells  calculation,  so  that  for  a given  number 
of  these,  R^  is  likely  to  be  smallest  and  the  last  output  table  of  thickness 
distribution  is  likely  to  be  nearest  to  that  giving  the  desired  upper-surface 
pressure  distribution. 

3.4  Fourth  problem:  hybrid 

In  this  problem,  the  upper-surface  pressure  distribution  is  again  specified 
everywhere,  but  different  second  conditions  are  imposed  inboard  and  outboard  of  a 
spanwise  section  fairly  near  the  root,  q = q*  say.  For  q > n*  the  wing 
thickness  distribution  is  specified.  In  the  course  of  iterations  on  this  part 
of  the  problem,  a doublet  distribution  2(S,n)  is  calculated  (for  q > q*)  and 
repeatedly  adjusted  as  in  the  second  problem;  it  is  then  necessary  to  extra- 
polate l in  some  way  to  the  inboard  region  0 < q < q*  , and  one  way  to  do 
this  is  simply  to  require  that  the  computed  inboard  doublet  distribution  be 
independent  of  spanwise  position: 

iU, q)  = n(5,q*)  (0  < q < q*)  . 

In  linear  theory  this  would  be  equivalent  to  maintaining  the  chordwise  load 
distribution  right  into  the  root,  and  in  our  problem  only  small  departures  from 
this  condition  should  result.  Finally,  to  satisfy  the  upper-surface  pressure 
condition  for  0 < q < q*  we  adjust  the  thickness  distribution  in  that  region 
as  in  the  third  problem. 
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The  program  has  been  arranged  to  treat  n = h*  as  the  third  collocation 
station  ourboard  from,  but  not  including,  the  root,  so  that  there  are  two 
collocation  stations  inboard  of  n = n*  on  which  the  thickness  distribution  is 
to  be  adjusted.  In  calculating  the  thickness  perturbations,  we  consider  that  it 
is  essential  to  ensure  that  the  thickness  distribution  for  q < n*  always  fairs 
smoothly  into  the  given  distribution  for  q > q*  , and  the  most  practical  way  to 
do  this  seems,  to  fit  (by  least  squares)  quadratic  curves  with  the  required 
(zero)  first  derivative  spanwise  to  the  calculated  thickness  perturbations  at 
each  value  of  £ . 

The  same  curve  fit  is  used  to  extrapolate  the  new  thickness  distribution 
to  the  root  q = 0 , which  is  again  not  a collocation  station  but  is  still  an 
output  station.  This  curve  fitting  means  that  near  the  root  the  upper-surface 
pressure  condition  is  now  satisfied  in  a mean  sense  only,  at  each  collocation 
station  value  of  £ . 

4 RESULTS 

4. 1 First  problem 

g 

Results  were  available  for  the  cambered  and  twisted  RAE  Wing  'B'  at 
Mach  number  0.8,  from  the  author's  direct  program'.  Wing  'B'  has  planform 
aspect  ratio  6,  taper  ratio  1/3,  straight  leading  and  trailing-edges  on  each 
half-wing,  and  mid-chord  sweep  angle  30°  (Fig. 2);  the  chordwise  thickness 
distribution  is  that  of  the  9 per  cent  thick  RAE  101  section.  As  a test  case 
for  the  first  problem,  the  final  output  planar  doublet  strength  from  these 
results  was  input  as  data;  its  behaviour  at  three  spanwise  stations,  near  the 
root,  in  mid-semispan  and  near  the  tip,  is  shown  in  Fig. 2a.  The  program  was  run 
for  four  iterations. 

As  the  planar  doublet  strength,  and  hence  the  velocity  field  ii  , is 
fixed,  we  may  expect  the  overall  spanwise  loading  properties  (which  depend  only 
on  the  doublet  distribution  in  linear  theory),  the  residual  error  R^  and  hence 
the  camber  and  twist  distributions,  to  settle  down  fairly  quickly  to  their  final 
values;  this  expectation  is  borne  out  by  Figs. 2b  and  2c  in  which  the  difference 
between  the  spanwise  twist  and  loading  distributions  calculated  at  the  first  and 
fourth  iterations  can  hardly  be  seen  on  the  graph.  The  corresponding  camber 
distributions  at  three  stations  are  shown  in  Fig. 3,  and  except  near  the  root, 
the  difference  between  results  for  the  first  and  fourth  iterations  is  also  very 
small.  Also  shown  in  Figs. 2b  and  3 are  the  actual  values  for  Wing  'B',  which  do 
exhibit  noticeable  differences  with  the  converged  results.  At  the  outboard 
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stations  in  mid-semispan  and  near  the  tip,  these  differences  can  be  attributed 
to  numerical  error  in  the  trapezoidal  integration  method  used  to  evaluate  the 
camber  and  twist  integrals  (3-3),  (3-4),  when  data  is  available  only  at  seven 
interior  Weber  points;  for  this  small  number  of  chordwise  points,  which  has  been 
taken  for  demonstration  purposes  only,  the  error  in  the  twist  distribution  is 
only  of  the  order  5 per  cent,  and  would  be  reduced  (and  more  detailed  input  and 
output  secured)  if  more  chordwise  points  were  taken.  The  large  difference  in 
camber  near  the  root  (in  Fig. 3)  is  almost  certainly  due  to  the  absence  from  the 
design  program  of  the  line  singularities  introduced  at  the  root  in  the  direct 
program;  although  these  line  singularities  have  their  major  effect  on  the  root 
section,  tney  do  have  some  effect  on  the  calculated  planar  doublet  strength 
(through  cross-coupling  with  the  boundary  conditions)  at  the  first  outboard 
station.  These  factors  lead  also  to  a considerable  difference  between  values 
from  direct  and  design  programs  in  the  upper-surface  pressure  C ^ near  the 
root,  as  shown  in  Fig. 4;  however,  for  the  two  outboard  stations,  notwithstanding 
the  slight  differences  in  camber  there,  the  differences  between  converged  values 
and  values  from  the  direct  program  could  hardly  be  distinguished  on  this  scale, 
and  have  been  omitted  for  clarity.  Fig. 4 shows  principally  that  C ^ takes 
longer  (but  not  unacceptably  longer)  to  settle  down  than  the  camber  and  twist 
near  the  root,  and  slightly  longer  at  the  outboard  stations;  this  is  because 
the  source  distribution  has  to  be  adjusted  repeatedly  to  reduce  the  boundary 
condition  error  R^  , and  we  know  from  experience  with  the  direct  program  that 
this  error  is  likely  to  be  largest  near  the  root,  and  has  a rather  uneven 
convergence  ratio. 

Results  from  the  program  for  this  first  design  problem  (to  which  we  shall 
refer  as  Option  1)  have  also  been  used  as  test  cases  for  the  other  problems,  a 
procedure  which  seems  likely  to  produce  consistent  results  as  the  basic  calcu- 
lation methods,  and  the  principal  source  of  error  in  these  demonstrative  cases 
(in  numerical  evaluation  of  the  camber  and  twist  integrals),  are  the  same. 

Also,  any  comparisons  with  results  from  the  direct  program  would  be  bedevilled 
by  the  effects  of  the  line  singularities,  just  as  we  have  already  seen  for  this 
first  problem;  by  instead  taking  results  from  the  design  program  as  test  data 
and  convergence  targets,  it  is  much  easier  to  assess  the  behaviour  of  the  other 
design  programs  near  the  root. 
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4. 2 Second  problem 

In  the  last  section  we  obtained 
0.8  from  the  design  program  Option  1. 
the  second  problem,  we  input  the  final 
distribution,  along  with  the  original 
program  for  four  iterations. 


results  for  RAE  Wing  'B'  at  Mach  number 
To  gain  some  experience  of  the  program  for 
upper-surface  pressures  C ^ as  a target 
thickness  distribution,  and  again  ran  the 


The  target  distributions,  and  the  results  from  the  first  and  fourth 
iterations,  are  shown  in  Figs. 5-7.  The  targets  are  the  same  as  the  final 
distributions  shown  in  Figs. 2-4,  and  are  now  represented  by  full  lines.  A 
first  glance  at  Figs. 5b  and  c for  the  twist  and  spanwise  loading  characteristics 
suggests  that  convergence  is  good  at  the  outboard  stations,  mid-semispan  and 
near-tip,  but  that  the  results  are  not  fully  converged  near  the  root;  Fig. 5a 
shows  that  the  doublet  strength  near  the  root  is  also  converging  rather  slowly. 
The  camber  distributions  (Fig. 6)  tell  a similar  story:  outboard,  very  good 
convergence  to  the  target  from  a rather  poor  first-iteration  result;  near  the 
root,  still  some  way  to  go  though  the  general  shape  of  the  curve,  including  the 
hump  near  the  root  trailing-edge,  is  well  predicted.  Fig. 7 shows  the  corres- 
ponding behaviour  of  C , and  it  is  rather  surprising  that  the  remaining 
change  required  in  C near  the  root,  in  particular  near  the  apex,  corresponds 
to  so  large  a remaining  change  in  camber  and  twist  according  to  .igs.5b  and  6. 
This  may  be  due,  amongst  other  things,  to  the  proximity  of  the  root-line 
singularity  in  upwash  corresponding  to  the  kinked  doublet  distribution,  so  that 
a small  change  in  doublet  strength  may  produce  a large  change  in  upwash,  a-d 
hence  in  camber  and  twist. 


It  is  perhaps  worth  commenting  on  the  first  guess  for  the  doublet  strength 

as  shown  in  Fig. 5a.  This  first  guess  is  calculated  using  the  older  version  of 

Lock's  method^  in  which  the  velocity  components  due  to  thickness  are  first 

£ 

estimated  using  the  RAE  Standard  Method  ; but  it  is  difficult  to  see  how  to 
improve  it  substantially,  as  the  accompanying  first  guess  at  the  source  distribu- 
tion is  such  that  the  first  estimate  for  C is  not  too  far  out  near  the  root, 

pu 

and  indeed  seems  excellent  towards  the  trailing-edge  for  this  case.  These 

results,  with  the  attendant  errors  in  the  first  estimates  of  camber  and  twist 

shown  in  Figs. 5b  and  6,  suggest  that,  in  the  design  problem,  second-order 

effects  are  very  important  near  the  root,  and  that  even  though  the  pressure 

coefficient  C from  a first-order  scheme  may  be  near  to  the  required  value, 

Pu 

it  is  necessary  to  check  that  the  boundary  conditions  on  the  wing  surface  are 
well  satisfied  by  the  velocity  fields  assumed  or  implied  by  such  a first-order 
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4.3  Third  problem 

It  may  happen  that  a wing  can  be  designed  to  support  a uniform  spanwise 
distribution  of  upper-surface  pressure  and  of  loading  near  the  root,  with  bene- 
ficial aerodynamic  consequences,  by  increasing  the  section  thickness/chord 
ratio  near  the  root.  If  this  design  also  avoids  an  unrealistic  increase  in  root 
twist,  and  happens  to  be  structurally  sensible  and  convenient,  so  much  the 
better. 

To  provide  a test  case  with  a known  solution,  a wing  somewhat  similar  to 
RAE  Wing  ' B'  and  designated  Wing  ' B',  was  designed  first,  using  Option  1.  The 
thickness  distribution  was  specified  as  the  RAE  101  section,  the  thickness /chord 
ratio  t being  9 per  cent  at  and  outboard  of  the  collocation  station 
1*  = 0.1563  , rising  parabolically  to  13.5  per  cent  at  the  root  (see  also 
Fig. 8a) : 
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V 0. 1 563/ 


(n  < 0.  1 563)  . 


The  behaviour  of  the  specified  doublet  distribution  at  three  stations,  including 
the  inboard  region  ri  < 0.1563  , is  shown  in  Fig. 8b.  (This  choice  of  doublet 
distribution  was  a historical  accident  based  on  earlier  work  on  the  fourth 
problem. ) 

The  results  from  this  run,  representing  the  target  distributions,  are 
again  shown  as  the  continuous  lines  in  the  remaining  figures.  The  program  was 
run  for  five  iterations.  Fig. 8c  shows  the  convergence  of  the  twist  distributions, 
and  we  see  that  the  target  is  nearly  attained  everywhere,  and  that  the  first 
shot  was  not  far  wide  of  the  mark  either.  (In  this  case,  the  effect  of  root 
thickening  on  the  design  root  twist  is  marginal:  a degree  or  so  less  than  the 
values  shown  for  Wing  'B'  in  Fig. 2b,  for  the  same  lift  coefficient  •) 

Fig. 9 shows  the  convergence  of  the  thickness  distribution.  Considering 
first  the  two  outboard  stations,  we  see  that  the  results  for  the  third  and  fifth 
iterations  are  virtually  identical,  and  that  the  corresponding  pressure  distri- 
butions (Fig. 10)  and  camber  distributions  (Fig. 11)  are  very  nearly  on  target, 
being  indistinguishable  except  near  the  leading  edge.  However,  the  converged 
thickness  distributions  are  not  quite  on  target.  This  must  be  due  to  numerical 
error  in  integrating  3Azt/3x  , like  the  corresponding  numerical  error  demon- 
strated for  the  camber  in  Fig. 3.  The  relative  error  seems  smaller  than  that  in 
the  camber;  a likely  mitigating  cause  is  that  part  of  the  error  in  z is 
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picked  up  when  the  velocity  components  are  evaluated  on  the  incorrect  surface. 
This  error  would  likewise  decrease  as  the  number  of  collocation  points  chordwise 
is  increased. 

Near  the  root,  the  thickness  distribution  converges  very  slowly,  as  the 
gap  between  third  and  fifth  iterations  in  Fig. 9 shows.  The  difference  between 
the  results  for  the  fifth  iteration  and  the  target  distribution  is  larger  than 
either  the  remaining  convergence  leeway  we  might  expect  or  the  numerical  error 
in  integration  we  could  anticipate  from  the  results  outboard,  but  can  easily  be 
imagined  as  the  sum  of  these  two  contributions. 

The  pressure  distribution  in  this  region  (Fig. 10)  also  exhibits  slow 

convergence  in  the  first  half-chord,  but  the  final  shortfall  of  the  target  is 

little  different  from  that  outboard.  It  may  be  an  inherent  difficulty  for  the 

third  problem  that  the  pressure  distribution  is  less  sensitive  to  the  thickness 

distribution  near  the  root  than  outboard.  It  is  not  easy  to  decide  how  to  cope 

with  this  difficulty.  An  over-relaxation  factor  of  about  2 could  be  introduced 

near  the  root  (it  is  not  needed  outboard),  to  speed  up  the  convergence,  though 

we  would  prefer  to  build  up  more  experience  before  citing  this  as  the  universal 

panacea.  It  would  also  help  if  the  basic  estimate  z _ could  be  improved  near 

tB 

the  root;  inspection  of  the  detailed  computer  output  shows  that  the  improvement 

obtained  with  Lock's  formula  is  not  remarkable,  so  that  the  improvement  would 

need  to  be  fairly  drastic;  we  also  note  that  the  basic  estimate  definitely 

overpredicts  the  thickness  at  the  mid-semispan  station,  whereas  the  thickness 

is  underpredicted  at  the  root.  (There  was  no  zero  in  the  first  estimate  for 

z , and  so  this  estimate  did  not  have  to  be  modified  as  described  in 
t B 

section  3.3  and  Appendix  C.)  There  is  also  the  possibility  of  setting  up  an 
inner  iteration  cycle  for  the  successive  increments  Az^.  , as  was  done  for  the 
doublet  strength  in  the  second  problem,  but  it  is  very  likely  that  in  this 
problem  the  dominant  disturbance,  not  taken  into  account  by  the  formula  (3-7), 
is  the  change  in  the  velocity  fields  of  the  two  singularity  distributions  when 
computed  on  different  thickness  surfaces,  rather  than  just  the  change  in  v 
and  due  to  perturbation  sources.  These  changes  could  be  estimated,  using 

Taylor  series,  if  the  program  were  rearranged  to  store  the  field  derivatives 
Su^/Sz  , etc.  near  the  root,  which  are  currently  overwritten  to  save  core  store. 

The  remaining  error  in  the  camber  distribution  near  the  root  (Fig. 11)  is 
obviously  associated  with  the  errors  in  the  other  field  quantities,  but  is  at 
least  an  order  of  magnitude  smaller  than  the  remaining  error  in  the  thickness 
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distribution.  We  might  expect  this  since  the  camber  and  twist  distributions 
depend  largely  on  the  input  doublet  strength.  On  the  other  hand,  the  outboard 
results  for  the  first  iteration,  which  do  differ  somewhat  from  the  target 
distributions,  show  that  the  agreement  in  the  twist  distribution  at  the  same 
stage  (Fig. 8c)  is  probably  fortuitous  and  that  second-order  effects  are  important 
for  this  problem,  outboard  as  well  as  near  the  root. 

4.4  Fourth  problem 

For  this  final  problem,  the  same  test  case  was  chosen  as  that  for  the 
third  problem:  the  wing  we  have  denoted  as  RAE  Wing  The  section  separating 

regions  where  different  conditions  are  applied  was  the  section  where  the  inboard 
rise  in  thickness /chord  ratio  begins,  n = n*  = 0.1563  . For  this  problem,  the 
basic  estimate  for  inboard  thickness  has  not  been  programmed  as  in  the  third 
problem;  instead,  the  first  guess  was  simply  taken  to  be  the  same  as  the  fixed 
outboard  distribution,  the  9 per  cent  thick  RAE  101  section.  The  program  was 
again  run  for  five  iterations. 

The  target  distributions  are  again  shown  as  full  lines  in  Figs. 12-14.  We 
see  that  in  mid-semispan  and  near  the  tip,  all  quantities  converge  well,  as  they 
did  for  the  second  problem,  to  which  this  hybrid  problem  is  essentially  equival- 
ent outboard.  We  also  observe  the  same  poor  nature  of  the  first  guess,  which 
corresponds  essentially  to  established  first-order  techniques,  despite  the  fact 
that  at  mid-semispan  the  chordwise  pressure  distribution  is  not  too  far  wrong. 

Near  the  root,  as  usual,  convergence  is  slow  and  the  thickness  distribu- 
tion (Fig. 12a)  has  not  converged  to  graphical  accuracy.  But  we  are  somewhat 
nearer  the  target  than  we  were  in  the  third  problem,  and  indeed  the  major  part 
of  the  remaining  error  could  be  just  the  numerical  integration  error  in  Azt 
(Fig. 9).  Convergence  of  the  root  twist  (Fig. 12b)  and  camber  (Fig. 13)  is  not 
quite  as  good  as  in  the  third  problem,  but  much  better  than  in  the  second  prob- 
lem. Here  too,  improvement  on  first-order  results  is  noted.  The  graphs  of 
section  lift  and  centre  of  pressure  duly  exhibit  the  expected  spanwise  invariance 
near  the  root  (Fig. 12c).  The  upper-surface  pressure  distribution  has  converged 
to  about  the  same  level  of  accuracy  as  in  the  second  and  third  problems  (Fig. 14); 
the  results  near  the  root  have  been  plotted  for  the  second  iteration  rather  than 
the  first  (for  which  only  the  ad  hoc  first  shot  for  the  thickness  distribution 
was  available)  to  show  again  the  level  of  error  when  only  one  thickness  pertur- 
bation is  calculated. 
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Wing  'B'  is  a special  case,  tailored  for  this  problem  so  that  the  root 
thickness  distribution  is  precisely  of  the  quadratic  spanwise  form  which  the 
inboard  least  squares  fit  demands  so  that  no  difficulty  arises  with  this 
smoothing  artifice.  Moreover,  the  interaction  between  the  unknown  doublet 
strength  and  the  unknown  thickness  inboard  happens  to  be  favourable,  so  that 
overall  improvement  in  convergence  is  found  compared  with  both  the  second  and 
third  problems.  Nevertheless,  the  results  for  this  hybrid  and  somewhat 
pathological  problem  at  least  seem  promising. 

5 CONCLUSION 

In  this  Report  we  have  studied  four  wing  design  problems,  for  two  of  which 
the  solutions  have  already  been  considered  in  principle  by  Weber**.  The  imple- 
mentation of  the  first  problem  (given  thickness  and  doublet  strength,  or 

first-order  loading,  £)  was  straightforward  and  the  program  for  it  (Option  1) 
converged  rapidly. 

The  second  problem  (given  z and  the  upper-surface  pressure  distribution 

C ) and  the  third  problem  (given  £ and  C ) have  been  satisfactorily  resolved 
pu  v 6 pu 

in  the  outboard  wing  regions,  mid-semispan  (more  or  less  sheared-wing  station) 
and,  perhaps  surprisingly,  near  the  tip;  but  they  have  proved  far  less  tractable 
near  the  root  of  a swept  wing.  For  the  second  problem,  the  camber  and  twist  are 
very  sensitive  to  the  doublet  strength  at  the  root,  as  we  might  expect  from 
first-order  theory,  and  there  is  considerable  cross-coupling  between  £,  C 
and  the  residual  error  R^  in  the  symmetric  boundary  condition.  Convergence 
has  been  secured  for  a swept  wing  at  high  Mach  number,  but  not  for  another  in 
incompressible  flow,  for  which  it  seems  likely  that  the  sidewash  and  upwash 
velocity  components,  which  are  not  reduced  in  scale  by  the  Prandtl-Glauert 
factor  relative  to  the  streamwash,  make  it  almost  impossible  to  find  a suitable 
approximate  doublet  distribution  to  satisfy  the  upper-surface  pressure  condition 
near  the  root  trailing  edge  in  the  first  one  or  two  iterations,  so  that  an 
optimization  technique,  yet  to  be  devised,  is  required.  For  the  third  problem, 

C does  not  seem  to  be  very  sensitive  to  z^  at  the  root,  again  Rt  has  to 
be  allowed  to  settle  down  before  z is  adjusted,  and  convergence  is  slow. 

Also,  since  the  velocity  fields  due  to  the  complete  source  and  doublet  distribu- 
tions have  to  be  computed  on  each  new  thickness  surface,  the  program  runs  for 
rather  longer  than  Options  1 or  2.  But  the  program  (Option  3)  has  not  actually 
failed  to  converge  for  any  case  for  which  a solution  is  known. 
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The  fourth  problem  is  a hybrid  of  the  second  and  third  problems  in  which 
l is  to  be  determined  outboard,  and  the  thickness  distribution  inboard  of  a 
certain  section  near  the  root.  In  our  opinion  this  hybrid  stands  an  excellent 
chance  of  giving  the  wing  designer  the  solution  to  his  commonest  problem  out- 
board, while  avoiding  the  severe  practical  difficulties  encountered  near  the 
root  in  Option  2 and  instead  determining  the  required  increase  of  thickness 
there,  to  maintain  good  flow  quality  and  to  give  additional  structural  strength. 
The  results  for  one  particular  case  show  better  convergence  than  those  from 
either  Options  2 or  3,  and  while  the  dangers  of  arguing  from  one  case  are 
realized,  this  option  (Option  4)  seems  promising  and  worth  bringing  to  the 
attention  of  designers. 

It  may  be  asked  why  we  have  used  the  programs  to  obtain  results  for  wings 

at  the  high  subcritical  Mach  number  0.8,  since  it  is  known  that  the  first-order 

Prandtl-Glauert  rule  can  only  be  expected  to  give  good  results  for  low  sub- 

critical  Mach  numbers  for  which  the  flow  nowhere  approaches  sonic  speed.  One 

answer  is  that  a design  application  can  be  envisaged  at  high  subcritical  Mach 

numbers,  if  a shockfree  flow  is  sought  and  if  a shockfree  solution  of  good 

quality  sufficiently  close  to  the  design  condition  is  available  from  another 

method,  or  from  experiments,  which  might  even  be  for  a wing-body  or  wing-nacelle 

combination.  In  this  case,  we  would  assume  that  changes  in  wing-body  interaction 

effects  due  to  small  perturbations  on  the  given  wing  are  negligible  over  the 

major  part  of  the  wing.  Let  us  denote  the  upper-surface  pressure  distribution 

from  this  given  solution  by  C _ and  that  required  by  C _ . Using  our 
, pu,E  M 7 pu,D 

direct  program  , we  can  compute  the  Prandtl-Glauert  solution  p for  the 

given  isolated  wing,  and  it  will  exhibit  an  error  C _ - C _ which  would 

pu,P  pu,E 

also  include  the  interaction  effects  in  a wing-body  combination.  Since  ^ 

does  not  differ  much  from  C „ , we  could  expect  that  the  wing  we  seek  does 

pu,E_ 

not  differ  much  from  the  datum  wing,  and  that  the  error  in  the  Prandtl-Glauert 

solution  Cpu  for  the  wing  sought  would  not  differ  much  from  the  error  in  the 

Prandtl-Glauert  solution  C _ for  the  datum  wing: 

pu,P 

C -C  n = C - C 
pu  pu,D  pu,P  pu,E 

hence 


= C _ + C - C 

pu,D  pu,P  pu,E 


I 


31 


If,  then,  we  design  our  wing  to  have  the  upper-surface  pressure  distribution 
, it  seems  likely  that  the  true  result  will  be  close  to  the  target  distri- 
bution C _ . Which  option  is  used  would  depend  on  whether  the  thickness  or 
.pu,D. 

the  loading  distribution  is  to  be  retained;  we  remark  that  if  the  thickness 
distribution  is  to  be  retained  outboard,  Option  4 might  be  better  than  Option  2 
(since  the  centre  line  distributions  in  the  direct  program  are  not  reproduced  or 
used  in  the  design  programs),  but  would  be  less  important  if  the  wing  thickness 
can  be  increased  at  the  root. 


It  is  also  worth  mentioning  briefly  a possible  use  in  landing  and  take-off 
design  studies.  The  programs  cannot  be  used  to  design  wings  with  separate  slats 
or  other  high-lift  devices,  but  they  might  be  used  to  obtain  a suitable  camber 
line  for  the  RAE  variable  aerofoil  mechanism  (RAEVAM) , which  can  under  some 
circumstances  compete  with  separate  high-lift  devices. 
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LOCK'S  FORMULA  FOR  THE  FIRST  ESTIMATE  OF  WING  LOAD  DISTRIBUTION 

Lock3  has  proposed  a formula  for  the  total  velocity  Q on  a wing  in 
compressible  flow,  suitable  for  use  in  the  design  problem: 


DQ 


8n  f / ,(3)  > 

COSa+F  UtB  COS  a 4 U*B  1 + V"  sec  Am 

n \ n ; 


v cos  a ± v „ 
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tB  JIB 
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0 “ - ^ sin2Ax  cos^a 


(A- 1 ) 


where  the  upper  signs  are  taken  on  the  upper  surface; 

a is  the  local  section  incidence,  and  is  therefore  equal  to  our  ; 
A^  is  the  local  sweep  angle  on  physical  wing; 

is  sweep  angle  of  maximum  thickness  line  on  physical  wing; 
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should  be  obtained  by  writing  B = 1 in  (A- 1 ) , giving  Q = "Q."  , and 
taking  "C  ."  = 1 - "Q."2  ; n 1 

6 pi  xi  ’ 

D - 1 + sec2A*  (pr)/&2n  * 

UtB(x»y.O),  vtB(x,y,0)  are  estimates  from  linear  theory  for  the  velocity 
components  due  to  wing  thickness,  for  example,  see  Ref.l; 

UJ,B^X,y’^’  v?B<x.y.°>  represent  the  thin-wing  components  due  to  doublets  in 
linear  theory,  and  are  to  be  determined; 


.(3) 


represents  the  second-order  interaction  effect  between  wing  thickness 
and  incidence. 
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It  is  convenient  to  mention  two  details  here.  First,  Lock's  formula  is  a 

semi-empirical  extension  of  the  Goethert  rule  to  take  into  account  higher-order 

compressibility  corrections  through  the  factor  B . Since  our  method  is  based 

n 

on  the  affine  transformation  which  yields  the  Goethert  rule,  there  is  no  point 
in  including  these  higher-order  terms  O^M^)  in  B , and  so  we  simply  take 
B^  = Sn  . This  reduces  (A-l)  to  equation  (3-5)  of  the  main  text. 

(3) 

Secondly,  we  have  slightly  modified  the  classical  derivation  of  S to 
derive  an  expression  which  seems  more  consistent  with  the  problem  at  hand  than 
the  one  given  in  Ref. 6,  even  though  it  is  only  one  of  the  several  second-order 
effects  present,  and  only  holds  for  two-dimensional  swept  wings.  Consider  such 
a wing,  with  uniform  sweep  angle  Ax  = A , and  incidence  a . The  analogous  wing 
will  have  the  uniform  sweep  angle  A where 

tan  A = B tan  A 


We  have,  for  the  thin-wing  velocity  components. 


JIB 


u„„  tan  A 

eB 


(A- 2) 


and  hence,  with  u„_  = u /8  , we  have  the  corresponding  relation  for  the 
X.B  xB 

analogous  wing  in  affine  space: 


VHB  = “ t3n  A 


The  first-order  boundary  condition  w = - a is  satisfied  by  the  doublet 

xB 

distribution  £_  with 
d 


= 1 


£B  e \ £ 

The  second-order  boundary  condition  gives: 


*b  = a(J~FS)"  cos  A ’ 


A hi  Ift  + lit  I^B 

.2  355  V£B  3y  Zt  3z 

P 


Using  the  sheared-wing  relation  3/3y  = - tan  A 3/3x  , and  the  zero  divergence 
of  in  our  affine  space,  we  obtain 
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cos  A sec^A 
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1 + sin  A ) W - T 25(1  - 5) 


c is  the  analogous  wing  chord.  Transforming  the  square  bracket  to  physical 
variables,  and  noting  that 


cos  A cos  A 
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sin  A 


sin  A 
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we  get 
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We  satisfy  this  boundary  condition  by  a doublet  strength  AS,  giving  the  further 
s treamwash 


1 


A0£  = IA*  = 
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(3) 

This  integral  differs  from  the  standard  definition  of  S by  the  presence  of 
8^  in  the  integrand. 

Using  (A-3)  again,  we  have  finally 


USB  + AUH  " U*B 


l + s(3*) 


and  similarly  for  v^  ; and  these  are  the  formulae  customarily  modified  to 
give  (A-l). 
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To  obtain  a first  estimate  for  u , Lock^  now  replaces  cos  a by  1 , 

X.D  T 

while  aT  is  still  an  unknown  quantity,  and  z^  by  z^  in  the  expression  for 
D ; and  the  dependence  of  v on  u is  taken  to  be  similar  to  that  exhibited 

XD  XD 

for  the  sheared  wing  by  (A-2): 


V*B  “ ' UJ,B  tan  Av 


where  tan  A = 
v 


2 

, cos  (AA  ) 

e2 5-2-  - * 

n cos  A 


= B 


cos  (AA  ) 
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2r 

cos  A 


(A- 4) 


Finally,  Q is  put  equal  to  the  design  value  Q . Then  (A— 1 ) becomes  a 

quadratic  in  the  variable  X = u _ 1 + S^  sec  A*/8  : 

S-B  ^ m n I 

AX2  + 2BX  + C = 0 


2 

where  A = sec  A 

v 

B = 1 + u „ - v „ tan  A 
tB  tB  v 

c = 0 + utB>2  + v2b  + (l  - *2)^°  " 1 ^ sin2Ax  - DQ  • 

Since  B is  of  the  order  (1  + small  quantities),  the  required  solution  is 


X = [(B2  - AC)^  - B]/A  . 


Hence  follow  the  estimates  for 


u^B  and  the  affine  doublet  strength: 


= 4u 


ZB 


= 4Bu 


ZB 


Since  a small  change  in  v has  a second-order  effect  on  the  value  of 

X D 

(A- 1 ) compared  with  that  of  a small  change  in  , we  can  discard  the  assump- 

tion (A-4)  and  replace  v^  in  (A-l)  by  the  value  calculated  from  the  first 
estimate  u and  the  more  accurate  formula': 

XD 


ZB 


_3_ 

9n 


c(n) 


! 


n)dC' 


\b  tan  A 


From  Z ^ we  can  also  compute  the  first  estimate  for  the  twist  and  write 


cos  aTB  in  (A-l).  This  leads  to  a revised  estimate  given  by 


cos  a 
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Lock's  formula  can  also  be  used  to  modify  an  initial  crude  guess  for  the 
thickness  distribution  z in  the  third  problem,  when  the  doublet  strength  is 

known.  From  the  doublet  strength,  we  can  find  u and  v . We  now  assume 

£B  36  B 

that  U£g  the  second  term  dominates  equation  (A- 1 ) , just  as  we  assumed  that 

UVB  same  term  dominated  the  equation  when  the  doublet  strength  was 

K • (3) 

unknown;  we  find  Vj-g»  ® and  S from  the  initial  guess  for  the  z 
distribution,  and  then  a revised  estimate  u^B  is  given  by 
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Appendix  B 

MODIFICATION  OF  PERTURBATION  QUANTITIES  WITH  MACLAURIN  SERIES 

The  main  iteration  cycle  involves  the  execution  of  the  lengthy  Ledger-Sells 
subroutine  to  obtain  the  velocity  fields  Aijt , Au^,  (or  ) accurately.  We 

might  obtain  a considerable  saving  in  computing  time  if  we  could  in  effect 
perform  one  iteration  cycle  by  using  suitable  estimates  for  these  velocity 
fields  instead;  we  would  then  hope  that  any  errors  due  to  the  estimates  would 
be  small  enough  to  be  absorbed  in  the  main  iteration  cycle  the  next  time  the 
velocity  fields  are  accurately  calculated. 

For  convenience,  we  drop  the  bracketed  iteration  superscripts  from  all 
quantities  except  the  residual  errors  R(=  ± Rfc  + R^)  . Let  us  suppose  that  in 
the  nth  iteration  cycle  we  have  computed  u = ii  + ti„  and  the  camber  and  twist 
distributions  z , , and  have  determined  perturbation  source  (and  possibly 

doublet)  distributions  Aq , AJt  and  perturbation  camber  and  twist  distributions 
Azg,  Aa^,  to  cancel  the  residual  field 


on  the  wing  surface 


z = z 


The  next  residual  field  on  the  perturbed  surface  z 
(suppressing  the  x,y-dependence) 


z + Az  will  be 
w w 


In  this  expression,  ii,  zw  and  are  first-order  small  quantities  while 

AG,  Az  and  Act_  are  (at  least)  second-order  small  quantities,  being  derived 
— w 1 . . 

from  such  expressions  as  R^  ' which  is  itself  second-order.  We  therefore 

expand  R^n’*^  in  powers  of  Azw  and  retain  only  quantities  up  to  and  including 

third-order,  ignoring  for  instance  terms  0(uz  Az  ) , 0(AuAz  ) . With  all 

w w w 

quantities  again  evaluated  at  z = z^  , this  leads  to 
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, (n , 1 ) 
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We  now  further  expand  Ci  in  powers  of  z (Maclaurin  series),  so  that  from 
this  point  on  all  quantities  are  evaluated  at  z = 0 , and  continue  to  retain 
terms  up  to  third  order.  We  can  ignore  further  contributions  from  the  second 
line  which  is  already  third-order.  Thus  only  Aw  contributes  further  to  the 
expansion: 


Aw(z  ) = Aw(0)  + z + 0/"z2Aw'\ 

w w 3z  \ w ) 

We  have  already  taken  our  perturbation  source,  camber  and  twist  distributions 
according  to  equations  (3-1),  (3-2)  to  make 


, \ , 3 Az 

» - a"<°>  *»ir 
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We  also  have 
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and  a similar  relation  for  Aw  . Hence  we  can  eliminate  w,  Aw 
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and  obtain 


(z  Av  + vAz  ) . 
w w 


Proceeding  on  the  lines  of  section  2,  we  write  for  upper  and  lower  wing  surfaces 


z=±z+z;  Az=±  Az^  + Az 

w t s w t s 
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with  similar  expressions  for  v,  Au,  Av  . 


problems,  Az  = 0 ; and 
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(For  the  first  and  second  design 
third  problems,  Au^  = Av^  = 0 .) 
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When  Azw  = 0 , these  expressions  reduce  to  those  given  in  Ref.l.  As  in  that 
document,  we  must  ensure  that  they  are  uniformly  valid  near  the  wing  leading 
edge  5=0.  We  require  at  worst 


Rt  = 0(5  5) 


E,  - 0(1) 


These  expressions  are  satisfactory  except  for  the  second  and  fourth  terms  in 

, (n, 1) 


the  square  brackets  in  R? 
this  square  bracket  (R*ot  say)  by 


3z 


We  introduce  a Riegels  type  factor  and  replace 


3Az 


Rc£  Aut  35  + °t  35  + 


\/('  * 6?) 


with 


3 z 3 Az 

\ = AQe  35T  + °£  TT  * 

This  completes  the  derivation  of  residual  errors  by  Maclaurin  series. 

To  estimate  these  residuals,  given  Aq  and  perhaps  AH  , we  need  quick 
estimates  for  four  velocity  components  on  the  chordal  plane  z = 0 . We  use  the 
approximations  for  A0t,  Av^  (derived  from  the  Standard  Method^)  and  AQ^,  Av^ 
as  set  out  in  equations  (5-12)  to  (5-17)  of  Ref.l.  We  also  need  , etc.  which 
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A 


are  just  the  accumulated  Au^  , etc.  on  the  plane  z = 0 . However,  instead  o 
simply  accumulating  them,  after  each  main  iteration  we  reset  ut  by  a Taylor 
series  based  on  the  accurate  value  at  z = z 

t ’ 


V°>  * Vzt>  - zt  JT 


and  similarly  for  v^,  , v^  . This  avoids  the  possibility  of  large  errors 

piling  up  in  the  accumulated  estimates  for  u^  , etc. 
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Appendix  C 

AN  IMPROVISED  FIRST  ESTIMATE  FOR  WING  THICKNESS  NEAR  THE  TRAILING  EDGE 

When  solving  the  third  design  problem,  if  the  first  estimate  for  the  wing 
thickness  takes  negative  values  near  the  trailing  edge  on  any  section,  we 
improvise  a distribution  to  replace  the  first  estimate  in  that  region 
temporarily,  until  the  main  iteration  scheme  gets  under  way.  In  such  a plane 
section  q = constant  , let  us  denote  the  first  chordwise  collocation  point 

A 

downstream  of  the  estimated  maximum  thickness  position  by  £ , and  denote  the 

• • A ' i > A A 

estimated  thickness  by  z and  its  derivative  by  v at  £ = £ . Then  we  would 
like  to  have  a curve  G(£)  which  passes  through  (£,z)  and  has  the  same  deriva- 

A A 

tive  G'  (£)  = v , becomes  rapidly  linear  further  downstream  and  vanishes  at 

A 

the  trailing  edge  £ = 1 . A function  which  becomes  rapidly  linear  for  £ > £ 
and  contains  three  unknowns,  a,  b,  c with  which  to  satisfy  the  other  three 
conditions  is 


G = ae_k(C  + b£  + c 


(£  < 5 < 0 


where  we  arbitrarily  choose  k so  that 


X = k ( 1 - £)  = 2 

so  that  the  exponential  term  dies  out  near  £ = 1 . X is  an  adjustable  program 
constant.  The  other  three  conditions,  in  order,  now  give: 

A A 

a + b£  + c = z 

A 

- ka  + b = v 
-X 

ae  + o + c = 0 

The  solution  of  these  equations  is  written: 

y(i  -0  + 2 

. -x 

1 - e - X 
b = ka  + v 
c = z - a - b£ 


We  observe,  as  a check,  that  a vanishes  and  G becomes  precisely  linear  if 

v = - z/(l  - £) 

a A 

which  is  the  slope  of  the  straight  line  joining  the  two  end  points  (£,z)  and 

0,0). 


3|  (3| 
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e* 

l 

M 

oo 

q 


Q 

R. 


u,  v,  w 
U 

—CO 

u 

U,  V,  w 
x,  y,  z 


SYMBOLS 

local  chord 

pressure  coefficient 

doublet  function  l sin  <f> 

planar  doublet  strength  (loading) 

free  stream  Mach  number 

planar  source  strength 

local  speed 

residual  error  in  symmetric  boundary  condition 

residual  error  in  antisymmetric  boundary  condition 

semispan  (taken  as  1) 

general  perturbation  velocity  vector 

general  perturbation  velocity  in  affine  space 

components  of  u 

free  stream  velocity  vector 

complete  velocity  vector:  + u.^  + u^ 

components  of  U 

local  Cartesian  coordinates  for  section 

leading-edge  ordinate 

wing  thickness  ordinate 

wing  camber  ordinate 

wing  section  ordinate:  z ± z 
0 st 


a 

B 

n 

ii 

n* 

A 

e 

a 


* 

a 


v 


local  section  twist 

wing  incidence  (taken  as  0) 

Prandtl-Glauert  factor:  iT]  - M^ 

* OO 

y/s 

y ; but  3/3n  denotes  differentiation  along  lines  of  constant  £ 

n = n*  is  section  dividing  root  and  outboard  regions  in  fourth 
(hybrid)  problem 

local  sweep  in  affine  space  (i.e.  on  analogous  wing) 
section  percentage-chord:  x = x^(y)  + c(y)£ 


angular  chordwise  coordinate:  E,  - JO  ~ cos  <}>) 

under-relaxation  factor  applied  to  estimate  of  output  by  vortex 
lattice  technique 


: 


Suffices 

B 

£ 

t 

u 


SYMBOLS  (concluded) 

basic  estimate 
due  to  doublets 
due  to  sources 
value  on  upper  surface 


Oversymbols 


design  quantities  (e.g.  C , Q ) 

pu  u 

quantities  in  affine  space  (e.g.  x = x/B,  u = uB) 
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Fig.  14  Iterates  for  upper-surface  pressure  distribution  on 

RAE  wing  "B",  M«  « 0*8.  Problem  4 
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